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^ Abstract 

,— I Community assembly is studied using individual-based multispecies models. The models have stochastic population dynamics 
with mutation, migration, and extinction of species. Mutants appear as a result of mutation of the resident species, while migrants 
have no correlation with the resident species. It is found that the dynamics of community assembly with mutations are quite 
^•different from the case with migrations. In contrast to mutation models, which show intermittent dynamics of quasi-steady states 

i_l interrupted by sudden reorganizations of the community, migration models show smooth and gradual renewal of the community. 

^5 As a consequence, instead of the 1 // diversity fluctuations found for the mutation models, 1 // 2 , random-walk like fluctuations are 
observed for the migration models. In addition, a characteristic species-lifetime distribution is found: a power law that is cut off by 

i-H a "skewed" distribution in the long-lifetime regime. The latter has a longer tail than a simple exponential function, which indicates 
an age-dependent species-mortality function. Since this characteristic profile has been observed, both in fossil data and in several 
pther mathematical models, we conclude that it is a universal feature of macroevolution. 

Ph Keywords: community assembly, macroevolution, migration, mutation, coordinated stasis, species-lifetime distribution, density 
Q dependent selection 
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Understanding the dynamics of biological macroevolution 
has been one of the most challenging topics in evolutionary 
theory. As more reliable and comprehensive fossil data are 
accumul ated, their st a tistica l aspects have attr acted increasing 



seen as selection in a dynam i cally changing f itness landscape 
(iGavrilets and Gravnen 1 1997tlGavrilets[ 120041) . This picture is 
in contrast with neutral models, in which no strong interspecie s 



interactions are included (Hubbelll l200lHPigolotti et alll2005l) . 



interest (Alrovetal., 2001; Newman, 2001). Meanwhile 



empirical studies of currently existing ecological com- 
munities provide snapshots of the evolutionary dynamics. 
Such fossil data and contemporary field data have previ- 
ously tended to be discussed separately. It has recently 
been recognized that issues of ecologic al and evo l utional 
scales can in fact be str ongly linked ( IThompsorl [1998, 
FT9991 lYoshida et al.L 120031) . Consequently, several models 
have recently been proposed to bridge the gap between 
ecological and evolutio nal timesc a les. These include the 
tangle d- nature model ( Hall et all 2002 ; Christensen et al 



In this article, the relation between temporal patterns on evo- 
lutionary timescales and community dynamics on finer ecolog- 
ical scales are studied with individual-based models. Simple 
stochastic processes, such as random walks and branching pro- 
cesses, have often been used for the interpretation of fossil data. 
However it is still an open question which aspects of evolution 
and extinction dynamics can be interpreted by such simple pro- 
cesses. In communities with complex interactions, do species 
compositions change gradually or inte rmittently showing coor- 



dinated stasis (DiMichele et al., 2004)? What are the necessar 



condi t ions for the "Red-Queen" hypothesis Jvan Valenl 11973 



Raud 1975; Doran et al 



20061: iFinnegan et all l2008t iBento 



20025 Idi Coll obiano et all 120031) , simplified versions 



of that model ([R ikvold and Zia, 2003; Zia and Ri kvold , 



2004 ISevim and Rikvoldl 120051: kikvold and Seviml 12007 



2009) to be valid? Do more realistic population-dynamics mod- 



Rikvoldfl2007t iFilotas et al.L I2010allbj) 
dCaldarelli et all Il998t 




the Webworld model 
2004 iMcKane . 



2001 



., Drossel et aL 
and others dChowdhurv et all 12003 : Shimada et al 
Tokita and Yasutomil 120031; iRikvoldlbOOgl) . These 



models consist of coupled population dynamics for each 
species, complemented by rules for introducing new species 
to the community. Survivability of individuals (or species) are 
determined by the totality of extant species. In other words, 
species undergo density dependent selection, which can also be 



els yield large avalanches as predicted in simplisti c SQC mod- 
els (IBak and Sneppenll 19931 iNewman and Palmeiil2003h ? An- 
swering these questions by theoretical model development will 
contribute to the interpretation of the statistics of fossil records 
and will also give insight into the conservation of currently ex- 
isting communities. 

Deliberately simple population-dynamics models with 
species turnover should answer questions about universal 
features of community assembly. The models considered 
here are simplified versions of the tangled-nature models 



(Rik vold and Zial 12003 
20071) . The tangled-nature model is 



Rikvold and Sevim, 2007; Rikvold, 



a simple individual- 
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based model, originally introduced by Ha ll and co-workers 
dHall et all l2002t IChristensen et all I2OO2I) and later simpli- 
fied bv IRikvold and Zial d2003l) . The evolutionary dynamics 
of various versions and derivatives of the tangled-nature model 
have been studied extensively. These previous studies have re- 
vealed that the model communities evolve intermittently rather 
than gradually. The evolution process consists of quiet pe- 
riods during which the species composition of the commu- 
nity rerrminsjieaTly COTistant, reminiscent of coordinated sta- 



1996h . 



sis (IDiMichele et all 120041) . and brief, active periods during 



which drastic rearrangement of the community composition 
takes place. The intermittency gives rise to approximate 1 // 
fluctuations of the species diversity, approximate r 2 species- 
lifetime distributions, and power-law duration distributions for 
quasi-steady states (QSSs). The 1 // fluctuations are distinct 
from those of the naive random walk: fluctuations on short and 
long timescales are self-similar and strongly correlated. Inter- 
estingly, these statistic al features are ob s erved not only for mu - 
tualistic communities ( Rikvold and Zia , 20031: Rikvoldi 20071). 
but also for several p r edator -prey models (IRikvold and Sevim , 



2007t iRikvoldl 120071 l2009h . Even though the network struc- 



tures dev eloped in these two types of communities are quite 
different (IRikvoldl 120071) . these dynamical statistics are simi- 
lar. Furthermore, neither model shows significa nt dependence 
on the choice of parameter sets (IRikvoldl 120071) . This implies 
that the observed dynamical features are universal with respect 
to the form of population dynamics, and supports the validity 
of such simple population-dynamics models. We expect that 
the dynamics on evolutionary timescales are characterized by 
a few key factors, and here we are aiming to identify some of 
these. 

In this article, we show that the way in which new species 
are introduced into the community plays a significant role for 
the stability of QSSs by comparing two models, called "mu- 
tation" and "migration." In the "migration" model, we con- 
sider community assembly via totally random migration of new 
species transported from another community. Exotic immi- 
grants are generally uncorrelated with native species since they 
have evolved in different environments. Therefore, any kind of 
new species may appear, regardless of the resident species. On 
the other hand, possible candidates for new species may often 
depend on the species composition of the resident community, 
for example in the case that the community is isolated and its 
evolutionary driving force is mainly sympatric speciation. This 
aspect is modeled in the "mutation" model. If the diversity of 
potential newcomers is limited, the community is expected to 
be more stable because it only has to be resilient against a lim- 
ited number of new species. This limitation is expected to have 
a significant effect on the stability of the community and its dy- 
namics on evolutionary timescales. In nature, some exotic inva- 
sions cause losses in the biological diversity of native species, 
and they are a major cause of de cline in global biodiversity 



dLodgelll993llFilotas et al.Ll2010allbh . ixotic invaders also have 



major economical consequences amounting to billions of dol- 
lars. Thus, this effect on ec osystem s is becoming an important 
focus of ecological s tudy (Lodgel 19931 Carlton and Gellei , 
19931 IVitousek et all 1 19971: iKolar and Lodgel l200ll ICarlton . 



2. Models 

Two forms of population dynamics (Model A and Model B) 
and two rules for the introductions of new species (migration 
and mutation) are investigated. In the following, we term the 
Model A (B) with the mutation (migration) rule "mutation (mi- 
gration) Model A (B)." 

2.1. Reproduction probability 

The models considered here are the simplifie d ver- 



sions of the tangled-nature model (|Rikvold and Zial [2003 
Rikvold and Sevimll2007l IRikvoldl 120071) . In these models, the 
population evolves stochastically in discrete, non-overlapping 
generations. Each individual of species / gives rise to F off- 
spring with a reproduction probability P[ before it dies. Oth- 
erwise it dies without offspring. The reproduction probability 
Pi for an individual of species / in generation t depends on the 
individual's ability to utilize the amount R of available external 
resources, and on its interactions with the population sizes n j(t) 
of all the species present in the community at that time. The 
form of Pi is taken as 



P,(R, {nj{t)}) = 



1 



l+exp[-A/(/?,{ny(r)})] ' 



(1) 



where 



^(R,{nj(t)}) = -bi + 



AW0 



y M u nj(t) _ Mot (2) 
4^ AWf) No' 



Here bj is the cost of reproduction for species / (always posi- 
tive), and 77/ is the ability of individuals of species / to utilize 
the external resource R. The matrix M defines the interactions 
between species. The total population size is N tM (t) = Yjj n AO, 
and No is an environmental carrying capacity that prevents 
Nt t(t) from diverging to infinity. The reproduction probabil- 
ity Pi(R, \n j(t)}) is a monotonically increasing function of A/, 
ranging over (0, 1). For a large positive A/ (small birth cost, 
strong coupling to the external resource, and more prey than 
predators), Pj approaches unity and the population of species 
/ increases. In the opposite limit of large negative A/ (large 
birth cost, weak or no coupling to the external resources, and/or 
more predators than prey), Pj goes to zero and the population 
decreases rapidly. We also note that the results shown below do 
not depend on the assumption that the number of offspring per 
individual is always F. For a reaso nable probabil i ty dist ribution 



of F, the results should be similar ( Murase et al. , 2010t) 



Two types of reproduction probabilities are considered in 
this article: Model A and Model B. Model A has no restric- 
tion on the form of the interaction matrix M. Therefore, each 
species makes various types of interactions with others, includ- 
ing predator-prey, mutualistic, and competitive interactions. 
Model B focuses on the energy transport through a food web, 
so the off-diagonal part of M is limited to being antisymmetric 
(M u = -Mji). Thus, if Mij > and Mj] < 0, then species / is 
the predator and J the prey. 



2 



In Model A, the reproduction cost bj and the external re- 
source R are zero; thus the first and the second terms of Eq. (f2]i 
disappear. The total population size is limited by the last term, 
which includes the carrying capacity No. The off-diagonal ele- 
ments of the interaction matrix Mjj are randomly drawn from a 
uniform distribution over [-1, +1], while the diagonal elements 
are set to zero. For Model A, F = 4 and No = 2000 are used in 
this article. This particular value of F is chosen such that per- 
turbations away from the fixed point in the single-species limit 
with vanishing mutation rate will dec ay monotonically , with- 
out oscillations or chaotic behavior (Rikvold and Zia, 2003; 
Filotas et all l2010albh As shown in dRikvold and Zial 120031 



Zia and Rikvoldi 12004 iRikvoldl |2007t iFilotas et all bOlOalbl) 



communities tend to evolve toward mutualism in Model A. 

In Model B, the external resource R is introduced. All the 
species have positive values of the birth cost bj, which are ran- 
domly drawn from (0, 1], and a certain proportion (0.05 is used 
in this article) of species can feed on the resource, i.e., the re- 
source couplings r]i are positive for primary producers or au- 
totrophs, and zero for consumers or heterotrophs. The resource 
R remains constant (here, 2000). The off-diagonal part of the 
interaction matrix is antisymmetric. Non-zero elements are as- 
signed randomly to Mjj (= -Mji) for I < J with probability 
c — 0.1, which is consistent wi th the connectance of food webs 
in nature dDunne et al. ■ l2002allbh . The nonzero elements of the 
interaction matrix are randomly chosen from a triangular dis- 
tribution on [-1, +1]. The diagonal elements of M, which rep- 
resent the intraspecies interactions, are selected randomly from 
a uniform distribution on [-1,0) for all the species. The envi- 
ronmental carrying capacity term is not included in this model 
(No = oo). The birth-cost term and the negative diagonal ele- 
ments Ma prevent species populations from growing to infinity. 
The fecundity F for Model B is set to 2 for the same reasons as 
discussed for Model A above. These models have fixed-point 
populations which can be cal culated exactly in the ab 



sence of mutat ions and migrations dRikvold and Seviml 12007 



Rikvoldi 12007). The dynamics of the populations are asymptot- 
ically stable but {«/} fluctuate around their fixed poin ts due to 
the demographic stochasticity (IZia and Rikvold , I2004 F1 

These models are advantageous since the fixed point and the 
linear stability can be analytically estimated. For the sake of 
comparison we used the same parameter sets as in earlier ar- 
ticles on the same models, but the results shown in the next 
section do not show qualitative differences for other reasonable 
parameter sets. 

2.2. Introduction of new species: migration versus mutation 

Introduction of new species and extinction of existing ones 
are essential for a macroevolution process of community as- 
sembly. New species are added to the system by "migration" 
or "mutation," which are different rules for the introduction of 



1 Strictly speaking, the unique stationary state for the mutation models is 
the state in which all species are extinct. However, the fluctuations around the 
locally stable fixed point, needed to cause complete extinction, are so extreme 
that the lifetime of a nonzero populat ion size should be Q(e N °) or 0(e R ), and 
thus for all practical purposes infinite (Zia and Rikvold. 20041) . 



new species. A species whose population becomes zero goes 
extinct. 

In the "migration" models, an individual of a new species 
appears in the community with properties (bj, r\j and Mjj) ran- 
domly assigned from the specified distributions. Thus, there is 
no correlation between the species existing in the system and 
the immigrants. Once a species goes extinct, it never appears 
again. The migrations happen at regular intervals of length t. 
This is definitely the simplest way to introduce new species into 
the system. 

In contrast to the "migration" model, in the "mutation" model 
the candidates for new species depend on the resident species. 
We consider the case that new species are limited to the ones 
that are phylogenetically close to the resident species. In order 
to include this effect, a coarse-grained "genome" space is intro- 
duced. Each individual has a bit-string "genome" of length L 
dHall et alll20q2tlChristensen et all \2002t. Idi Collobiano et aT 



20031: lGavriletsLl2004 . Each bit sequence corresponds to a dif- 



ferent species, which has a different phenotype from the neigh- 
boring species. Thus, the total number of potential species is 
2 L . All the values of bj, t]j, and Mjj are predetermined for each 
genotype at the beginning of the simulation. In every genera- 
tion, a mutation may happen to the genomes of the offspring: 
all the genome bits of appearing offspring, which amount to 
N ot xL, flip independently with a probability /u/L, so that the av- 
erage number of mutations per offspring individual is fi. These 
flips result in the appearance of new species and can be seen 
as a random walk along the ed ges of the hypercube defined 
by the states of the L-bit genome ( IGavrilets and Gravnerlll997l 



Gavrilets, 2004) 



The mutation rate per species, fi, determines how frequently 
new species appear. Thus, an individual moves randomly to 
a neighboring site in the L-dimensional hyper-cubic space by 
a mutation. The probability of m-bit mutations in a single in- 
dividual is small, 0(f/"). The properties of species / (bj, rjj, 
and Mjj) have no correlation with those of its neighbor species. 
This is a clearly idealized aspect of the model since phyloge- 
netically related species usually have phenotypic correlations. 
However, this model works as a good starting point for com- 
parison with the migration model because we can focus on the 
effects of the phylogenetic correlation. In particular we note 
that the dynamic insensitivity of mutation Model A to correla- 
tio n between the traits of pare nts and mutants was demonstrated 
bv lSevim and Rikvoldi J2005I) . 

Although the mutat ion model can be considered as a quasi - 
species model as well (lEigenll 1977tldi Collobiano et alll2003l) . 
the mutation used here does not necessarily refer to a point mu- 
tation of an actual genome when the model is discussed in the 
context of community assembly. The genome space is intro- 
duced in order to express the phylogenetic distance between 
two species and a bit-flip of the coarse-grained genome may 
correspond to a big leap in an actual genome space. The muta- 
tion model fits better than the migration model for cases where 
new species are always phylogenetically close to the resident 
species. Communities in isolated islands or lakes where sym- 
patric speciation is the major driving force of evolution may be 
relevant examples. On the other hand, communities geographi- 



3 



cally neighboring to a large species pool, such as communities 
in a continent, may be better described by the migration model. 
A detailed description of the speciation process is not included 
in either model. 

The crucial difference between the migration and mutation 
models is that the mutation model limits the variety of new 
species accessible from a given community. The reason why we 
use a binary string instead of a set of real numbers is that a bit- 
string maximizes the potential number of species while keeping 
the number of neighbor species reasonably small. Models with 
more alleles would be expected to show intermediate proper- 
ties between the results for the migration and mutation models. 
Although it is hard to estimate how many neighboring species 
should be assumed, the comparison of these two extreme mod- 
els highlights the effect of this limitation. Before adopting more 
complex models, we therefore start from these simple cases. 

3. Results 

3.1. Typical time series 

First we show typical time series of kinetic Monte Carlo sim- 
ulations for the four models. Genome length L = 25 and 22 are 
used for mutation Model A and B, respectively. Therefore, the 
numbers of potential species are 2 25 and 2 22 , respectively. The 
results shown in the following do not depend significantly on 
the precise values of L. The mutation rate p = 0.001 is used for 
the mutation models, and the migration interval r = 1 is used 
for the migration models. This mutation rate is selected for 
computational feasibility although it might be rather high com- 
pared with actual ecosystems. We confirmed that the power 
law behaviors shown below remain similar for mutation rates 
as small as 10~ 4 , and intermittency was qualitatively observed 
for p = 10~ 6 . The dynamic robustness of of mutation mod- 
els A and B with respect to the model parameter s (L, p., No, 
R) is discussed in Appendix C of iRikvoldl (12007b . In partic- 
ular, it is shown there that the diversities depend only sublin- 
early on Hubbell's fundamental biodiversity number, 2Rp or 
2Nop, which is proportional to the number of m utant organisms 
produced in each generation (IHubbelll 1200 lb . It is therefore 
justified to use large mutation rates in combination with small 
population sizes to overcome computational limitations. 

In this ar t icle, th e exponential Shannon-Wiener diversity in- 
dex (IKrebsl 1 1989b is used as a measure of biodiversity. This 
index is defined as the exponential function of the information 
entropy of the population distribution, D(t) = exp[S ({«/(f)!)L 
where S ({n 7 (t)}) = - 2{/|p ; ©>0} pi{t) In p/(f), with p/(f) = 
«/(f)/Mot(0- We adopted this measure in order to filter out un- 
successful mutants or migrants which have tiny populations and 
rapidly go extinct. 

Figure \T\ shows the dynamics of the diversity indices and the 
total population size for typical simulation runs for each of the 
four models. The mutation models show intermittent behav- 
iors, consisting of active and quiet periods. During the active 
periods, the diversity and the total population size show larger 
fluctuations, and the species composition changes quickly. Dur- 
ing the quiet periods, the species composition remains nearly 



constant, and the system is considered to be in a quasi-steady 
state (QSS). The community assembly proceeds intermittently, 
rather than gradually. 

In contrast, the intermittent behaviors are hardly observable 
in the migration models. It is difficult to find QSSs, at least 
from these figures. The system always fluctuates actively, indi- 
cating continuous renewal of the species composition. Thus, 
the assembly dynamics for the migration models are signifi- 
cantly different from the mutation models. A key ingredient 
of the intermittency in the evolution dynamics is the introduc- 
tion of the genome space, which limits the number of different 
mutants that can be produced by a given community. 
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Figure 1 : (Color online) Typical time series of exponential Shannon-Wiener 
diversity index and total population size for mutation and migration Models A 
and B. Mutation rate /i = 0.001 is used for the mutation models, and migration 
interval r = 1 is used for the migration models. Genome length L for mutation 
Model A and B are 25 and 22, respectively. 



3.2. Statistics of dynamical behaviors 

To evaluate the intermittency quantitatively, we calculated 
several statistics of the dynamical behaviors for each model. 
We performed simulations of 2 25 generations with 2 22 gener- 
ations as a "warm-up" period for Model A, and simulations 
of 2 26 generations with 2 24 generations warm-up for Model 
B. These warm-up periods are long enough to realize statisti- 
cally stationary states. For each model and species-introduction 
mechanism, the data were averaged over six independent runs. 
Each simulation run was started with a single, randomly cho- 
sen species (producer species for Model B) with a population 
size of 100 individuals. The details of this initial condition are 
totally insignificant, and the systems were completely "thermal- 
ized" during the initial warm-up periods. 

First, we estimated the QSS durations and calculated their 
probability density functions (pdf). One way to identify QSS is 
to introduce a cutoff on the logarithmic derivative of the diver- 
sity, dS jdt. The pdf's of this quantity (averaged over 16 gener- 
ations) are shown in Pig. El F° r all the models, each distribution 
has a shar p peak around the center and relat i vely w ide wings in 
both tails dRikvold and Zial 120031: IRikvoldl 12007b . The sharp 
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Figure 2: (Color online) Normalized histograms of the logarithmic derivative 
of the diversity index for Model A (a) and Model B (b). Mutation rate /j = 0.001 
is used for the mutation models, and migration interval r = 1 is used for the 
migration models. 



Figure 3: (Color online) Normalized histograms of the QSS durations for 
Model A (a) and Model B (b). Mutation rate // = 0.001 is used for the mu- 
tation models, and migration interval r = 1 is used for the migration models. 
The thresholds for estimating the QSS durations are 0.015 and 0.01 for Model 
A and Model B, respectively. 



peak around zero represents that the community is in a quiet 
period, while the large wings represent large rearrangements of 
the species composition. QSSs are estimated as the periods be- 
tween times when \dS /dt\ exceeds a cutoff. We adopted 0.015 
and 0.01 as the cutoff for Model A and Model B, respectively. 
We note that the migration models have larger wings than the 
corresponding mutation models. This indicates that the migra- 
tion models spend more time in active periods. 

The pdf's of the QSS durations are shown in Fig. [3] The 
QSS duration distributions for the mutation models show ap- 
proximate power laws. Approximate r 2 and f _I behaviors are 
observed up to 10 7 generations for mutation Model A and B, 
respectively. Thus, the community-assembly dynamics for the 
mutation models have long-time correlations. On the other 
hand, the migration models show faster decays and deviations 
from the power laws at certain characteristic time scales. For 
migration Model A and B, it starts to decay faster than a power 
law at about 10 3 and 10 4 generations, respectively. While the 
distributions are well approximated by power laws up to the 
characteristic time scales, another trend, which appears concave 
in a log-log plot, starts to emerge at longer time scales. Thus, 
the QSS distributions for the migration models have character- 



istic time scales, above which the distributions decay faster than 
power laws. This concave curve is not fitted well by a simple 
exponential function, neither for migration Model A nor B. We 
discuss a fitting function for this concave curve, known as a 
^-exponential, in the Appendix. 

The shorter QSSs for the migration models are also observed 
in power spectral densities (PSDs) of diversity indices. PSDs 
calculated for each model are shown in Fig. [4] For the mu- 
tation models, approximate l/f fluctuations (flicker noise) are 
observed over several decades. This indicates that the evolution 
dynamics for the mutation models have quite long time corre- 
lations. On the other hand, for the migration models, the PSDs 
saturate at low frequency. This indicates that the temporal be- 
havior of the diversity is uncorrelated for t > 1 /2nf c , where f c 
is the characteristic frequency at which the PSDs saturate. At 
/ > f c , the PSDs show approximate power laws 1 ff a . The ex- 
ponent a for the migration models are larger than for the muta- 
tion models, and are close to 2. Thus the dynamical behavior for 
the migr ation models are analo gous to an Ornstein-Uhlenbeck 
process (ICox and Millen, 1 19651) . in which the diversity index 
performs a random walk in an attractive "potential." The op- 



5 





10 11 F 




1 o • 






en 


10 9 ■ 


"O 




TO 


10 8 ; 


O 




pe 


10 7 - 


cn 


wer 


10 6 \ 


o 




Q- 


10 5 \ 




10 4 - 





■ 1 r- 


1 1 ■ — H 


i — ■ ■ — n ■ ■ — n — ■ ■ — n ■ ■ — r 

mutation ■ — i — ■ 
migration (x=1) — -x— ■ .... 
migration (t=64) -- »— ■ 
1/f 


H 










TK^ X, 1 

V *x 


1/ 


f 1.8 






.-*-x-x 


-x- x -x- 


*-x^x 




































(a)r. 


>/lodel / 


\ 








< 



10" 8 10" 7 10" 6 10" 5 10" 4 10" 3 10" 2 10" 1 
Frequency (1 /generation) 



10 u 



o 
cz 

Z3 



10"' 
10 



-4 



10 



-6 



CD 



10- 



& 10 



-Q 

CO 

O 



10 



n" 10 
10 



10 



12 



14 



16 



tt ■ ■ — n 


1 ■ — H 


— i ■ — n 


■ 1 — w 


t — • • — t — • ■ — n 

mutation 1 — 
migration 
1/t 2 


x- - 














































































(a) Model A 











10° 10 1 10 2 10 3 10 4 10 5 10 6 10 7 
Species lifetime (generations) 





10 13 r 




10 12 ■ 


ISU 


10 11 r 


CD 




"O 


10 10 f 


TO 


O 


10 9 ; 


CD 


o_ 




cn 


10 8 ' 


wei 


10 7 ■ 


o 


Q_ 


10 6 r 




10 5 ■ 



10 



mutation • — i — • 
migration (t=1) 
migration (t=1 6) -— x— • 
1/f 3 - 
1/f 1 




10" 7 10" 6 10" 5 10" 4 10" 3 10" 2 10" 1 
Frequency (1 /generation) 



o 
o 

Z3 



CD 



10 u 

10" 2 
10" 4 
10" 6 
10" 8 



10 

CO 

o 



10" 
10" 
10" 











m 


utation 1 — i — 1 
gration — x-- ■ 
1/t 2 - 3 










mi 
























^^x-. x 


x x 


















x 






















































ivioae 


I B 













10° 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 
Species lifetime (generations) 



Figure 4: (Color online ) Power spectral density of the diversity index for Model 
A (a) and Model B (b). Mutation rate fi = 0.001 is used for the mutation models, 
and several values of r are used for the migration models. 



Figure 5: (Color online) Species-lifetime distributions for Models A (a) and 
Models B (b). Mutation rate /i = 0.001 is used for the mutation models, and 
migration interval r = 1 is used for the migration models. 



timum of this potential corresponds to a balance between ex- 
tinctions and the invasion of new species. The characteristic 
frequency for migration models is approximately proportional 
to r -1 . The exponent a for small r is smaller than 2, but this 
may be due to a transient from l// 2 behavior to the saturation. 
We also calculated the PSDs of the total population sizes (not 
shown) and found the behaviors similar to the diversity indices. 

Figure|5]shows pdf 's of the species lifetime for the four mod- 
els. Both mutation Model A and B show approximate r 2 power 4 Discussion 
laws that continue up to around 10 7 generations. Thus, species 
that live quite long appear relatively frequently for the mutation 
models. However, the migration models show a different pro- 
file. Although the distributions show power-law decays at small 



t, they cross over to skewed profiles (Appendix Appendix A 1, 
which look similar to those we see in the QSS duration distribu- 
tions, i.e., concave on a log-log scale and convex on a semi-log 
scale. The shorter species lifetimes for the migration models are 
consistent with the shorter QSS durations since the rearrange- 
ment of communities cause large extinctions of species. 

We note clear similarities between the distributions of the 
QSS durations and of the species lifetimes for Model A. It im- 
plies that the lifetimes of long-living species are mainly deter- 
mined by the end of the QSS in which they live. In contrast, for 



Model B, the species-lifetime distributions are narrower than 
the QSS duration distributions. It means extinctions of long- 
living species frequently happen during each QSS. 

In conclusion, the migration models have a characteristic 
time scale in their dynamical behaviors, while the mutation 
models show power-law statistics. The differences between mu- 
tation and migration models are summarized in Table [1] 



Several studies of plant and animal assemblages from fossil 
records have revealed long-term persistence punctuated by pe- 
riods of rapid change (jPiMichele et al. , 2004 ; Brett and Bairdi 



1995; Brett et al., 1997). Such intermittent patterns of commu- 



nity assemblages, termed coordinated stasis, have been found 
on the basis of fossil records although counterexamples are also 
common. The possible origins of such intermittent patterns are 
still an open question. One of the aims of the present study is 
to find conditions under which coordinated stasis can arise in 
models of macroevolution. 

In previous studies of tangled-nature models with mutation 
dHall /et al.ll2002tlChristensen et all l2002t Idi Collobiano et al 



20031: iRikvold and Zial 120031: iRikvold and Sevimi 12007 
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Table 1 : Summary of the results for the four models studied. Representative behavior of QSS duration distributions (QSS), PSDs of the diversity and total population 
sizes (PSD), and species-lifetime distributions (SLD). 





mutation 


migration 


QSS 


1 It 1 (Model A), l/t (Model B) 


l/t 2 (l/t) + skewed profile 


PSD 


1// 


I jf 2 which saturates at low frequency 


SLD 


l/t 2 


l/t 2 + skewed profile 



Rikvoldl 120071) . it was found that the evolution proceeds 
intermittently rather than gradually. The systems spend most 
of their time in QSSs in which the species composition does 
not change significantly. The QSSs are interrupted by brief 
active periods, during which relatively large-scale community 
rearrangements happen. Duration distributions for QSS com- 
munities and species-lifetime distributions both show power 
laws. 

However, such QSSs become shorter and less well defined 
for migration models. The community-assembly dynamics is 
well described by an Ornstein-Uhlenbeck process, which is 
qualitatively different as it has a characteristic timescale. Thus 
it appears that the introduction of a genome space may play a 
significant role for the sustainment of QSSs. It is remarkable 
that mutation considerably differs from migration, even when 
the elements of the interaction matrix for the mutation models 
have no correlations with each other. The biggest difference 
between the mutation and migration models is the variety of 
possible new species. For the mutation model, the number of 
possible mutants is roughly proportional to L x D*, where D* 
is the number of major resident species. If most of the candi- 
dates are not suited to survive in the existing QSS community, 
the species composition does not change, and the emergence 
and rapid extinction of unsuccessful mutants is just repeated. 
In this case, a successful species can emerge by a two-bit mu- 
tation of a large-population species or by a single mutation of 
a low-population species, both of which happen with very low 
probability. Therefore, species which are neighbors in genome 
space of resident species act as a "protection zone." Intermit- 
tent dynamics in the community assembly can arise from the 
resulting limitation on the diversity of mutants. 

Intermittency has often been discuss ed for complex popula- 
tion dynamics of a fixed set of s pecies (jGavrilets and Hastings , 



1991' 



. Ives and Janseni . 
Ferriere and Cazellesl I1999Tk We note that the intermittency 



19981 iHuisman and Weissind, 12001 



observed in the present study is different from those studies 
since the population dynamics itself is asymptotically stable. 
The intermittency is observed at the level of the rearrangement 
dynamics of community compositions, not at the level of the 
population dynamics. 

We believe that the behaviors for the mutation models gradu- 
ally become similar to the migration models when the genome 
length L or the number of alleles increases. This is because the 
communities then are expected to have more paths to escape 
from the QSSs. For isolated environments such as lakes or is- 
lands located far from another community, the mutation models 
would fit better, while the migration models would fit better for 
communities adjacent to large species pools. 



Both for the migration and the mutation models, the species- 
lifetime distributions have heavier tails than simple exponen- 
tial functions. Thus, in communities with complex inter- 
species interactions, a naive Red-Queen hypothesis is not ex- 
pected to hold. The lifetime distribution estimated from the 
fossi l data also shows a heavier tail than a simp le exponen- 
tial (IFinnegan et all 120081 : IShimada et all 12.0031) . (See also 



Fig. IA.8l ) This fact indicates that extinctions due to the species 
interactions may explain the statistics in fossil records. 

The models we studied show a fundamental difference from 
the SOC models. Although large scale rearrangements of com- 
munities are observed (Fig. [2]), the extinction-size distributions 
do not show power laws. Therefore, these avalanches of extinc- 
tions are not critical processeses. Thus, we speculate that the 
key mechanism of the power laws observed for the mutation 
models are different from those of SOC models. 

In a limit that all interactions are equal to zero and all bj 
(and 77/) are equ al, the models become similar to neutral models 
(IRikvoldl, 120071) . In this limit, the species dynamics are charac- 
teristic of an off-critical branching process. We have prelimi- 
nary results for the neutral versions of the models. The species- 
lifetime distributions show a t~ 2 power law with an exponential 
cutoff, which agrees with t he return-time distrib ution of an off- 
critical branching process (IPigolotti et all 120051) . Intermittency 
is not observed for these models since no collective behavior 
between species is included. 

Since a large range of randomly distributed interactions and 
species-specific parameters (like tjj and bj) are available, the 
population dynamics selects those combinations that give rise 
to quasi-steady states (i.e., quasi-stable communities). In model 
A, this yields compact, mutualistic communities. The restric- 
tion to antisymmetric interactions in Model B limits the possi- 
bilities to predator-prey communities. The nonlinear character 
of the reproduction rate in our models make them more flexi- 
ble to portray a wide variety of different communities and their 
dynamics simply through the evolutionary selection of realized 
parameters, than is the case for traditional Lotka-Volterra mod- 
els. 

5. Summary 

Four types of biological community-assembly models were 
studied. By comparison between mutation and migration, we 
clarified that the limitations on the variety of possible mutants 
resulting from the introduction of a genome space plays a key 
role in the emergence of QSSs and intermittent evolution dy- 
namics that show power-law statistics over several decades. 
This also indicates that exotic migrants can destroy QSSs more 
easily than mutants of resident species. 
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The species-lifetime distributions for the migration models 
show a robust pattern: an approximate r 2 power law, consis- 
tent with a stochastic branching process, which crosses over 
to a skewed profile. This skewed profile can be reasonably 
fitted by a ^-exponential function, which is derived from an 
age-dependent species-mortality function. We believe that this 
skewed profile can be a natural candidate for the fitting of 
species-lifetime distributions in addition to simple exponentials 
and simple power-law functions. 

We emphasize the robustness of our results. Even though 
Model A and B have significantly different types of inter- 
actions and as a result develop different network structures 
dRikvoldl 120071) . significant similarities are found between the 
models. Furthermore, the profile of species-lifetime distribu- 
tions is similar to other models with quite different popula- 
tion dynamics (jChowdhurv et al. . 2003 : Shimada et al. . 2002 : 
Laird and JensenT 20061) . Further studies on both realistic and 
simplistic forms of population dynamics, such as the Webworld 
model, are planned. 
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Appendix A. Fitting with a ^-exponential function 

Since the migration models have characteristic time scales in 
their community-assembly dynamics, we expect that these are 
related to the migration interval r. Species-lifetime distribu- 
tions for migration Model A with several values of t are shown 
in Fig. |6(a)| The skewed profile in the longer- time regime shifts 
to the right as r increases, although the initial r 2 power law 
does not show significant dependence on r. The amount of the 
shift is approximately proportional to t. Thus, the distributions 
reasonably collapse onto a single curve for several values of t 
by the rescaling of time by t (Fig. |6(b)] i. Therefore, the typical 
time scale is determined by the number of migrations. 

The same plot for migration Model B is shown in Fig. |6(c)| 
The data scales reasonably with r although there are small de- 
viations. We speculate that these deviations are caused by the 
stochastic population fluctuations. 

QSS duration distributions normalized by r for migration 
Model A and B are shown in Fig. IA.7l The skewed profiles col- 
lapse well onto a single curve although the scaling is less clear 
for Model B. We also calculated PSDs for several values of t 
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Figure A. 6: (Color online) (a) Species-lifetime distributions P(t) for migration 
Model A with several values of r. (b) Species-lifetime distributions normalized 
by t, tPU/t), for migration Model A. It is fitted by a (/-exponential function 
with q = 1.296 and fi = 0.0208. (c) Species-lifetime distributions normalized 
by t, rP(r/r), for migration Model B, fitted by a (/-exponential function with 
q= 1.26 and /3 = 0.00018. 



and confirmed that the characteristic frequency f c where the sat- 
uration occurs is proportional to 1 jr. These are consistent with 
each other. Thus, it is confirmed that the time scale for the mi- 
gration models should be counted by the number of migrations. 
This is reasonable since the dynamical systems we studied are 
asymptotically stable and few major extinctions happen without 
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Figure A. 7: (Color online) QSS duration distributions normalized by t, 
P(//t)t, for (a) migration Model A and (b) migration Model B. A ^-exponential 
function with q = 1.334 and ji = 0.043 is included in (a) as a guide to the eye. 



being initiated by the addition of new species. 
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Figure A. 8: Lifetime distribution of families estimated from fossil records after 



IShimada etai]<2003l) . The original data consists of the list of presence of each 
families (present, absent, or unknown) in each geological era. Open circles and 
filled boxes represents the estimation by regarding the unknown case as present 
and absent, respectively. The line shows a ^-exponential fit with q = 1.23 and 
B = 0.03. 



2003): 



exponential function (IShimada et al. 

P(t) oc [1- (1 - q)pt\ m ~ q) , 



(A.l) 



where q and ft are fitting parameters. The asymptotic form of 
this ^-exponential function for t — > oo is a simple power-law 
decay, t l ^ l ^'\ In the limit of q — > 1, the ^-exponential function 
becomes a simple exponential function, e^'. The parameter 
j6 represents the inverse of the characteristic time scale. This 
function is derived from an age-dependent species-mortality 
function as shown below. We fitted a ^-exponential func- 
tion to the species-lifetime distribution for migration Model 
A with r = 256 for t > 2t, and obtained q = 1.296(26) 
and ft = 0.021(4) [number of migrations -1 ]. The agreement 
of the fitting func tions with the simula tion data is quite good. 



A previous study (IShimada et all 120031) revealed that the lifes 



pan of families in fossil data forms a similar "skewed pro- 
file," and the distribution is in reasonable agreement with a q- 
exponential function with q = 1.234(33) and /? = 0.0301(44) 
[million years -1 ]. (See Fig. IA.8H We note that q for this simu- 
lation model is in reasonable agreement with the q for the fos- 
sil data. We can compare the species-lifetime distribution ob- 
tained from the simulation with the fossil data although it may 
be debatable whether the comparison is valid because the fossil 
data show the lifespan not of species but of families. Com- 
parison of /? gives r = 0.7 million years for migration Model 
A, and r = 6 thousand years for migration Model B. From a 
fitting to the QSS duration distribution for migration Model A 
with t = 256 on t > 8t, q = 1.334(21) and p = 0.043(9) 
[number of migrations -1 ] are obtained. This q is also in rea- 
sonable agreement with the q for the fossil data. (The distri- 
bution for migration Model B is not clear enough for such fit- 
ting.) This pattern is so robust that it i s also expected to be ob 
served in other models, as indeed it is ( Chow dhurv et all 12003 



Shimada et al. , 2002 ; Laird and Jensen , 20061) 

The origin of the initial power laws in the species-lifetime 
distributions should be explained by the dynamics of unsuc- 
cessful species which fail to have a stable positive population. 
We expect that the population of unsuccessful species is de- 
scribed by a Galton- Watson like stochastic branching process 
since the population undergoes a reprod uctive process. The r 2 



decay is consistent with such a process dPigolotti et all 120051) . 

It is known that the <7-exponenti al function is obtained as a 
result of age-dependent mortality dShimada et all 120031) . By 
integrating Eq. ( IA.1I ). one can obtain the ratio of species that 
live longer than f, S (f), which is another ^-exponential function: 



S(t) 



r 



PQfW = \l-(l-q')0i\ 



1/(1-9') 



(A.2) 



where q' = 1/(2 - q) and P' = (2- q)p. The integral in Eq. ( lA2b 
converges only if 1 < q < 2, and this condition is satisfied for 
the simulation data. The probability that a species which has 
lived for t generations goes extinct in the (t + l)th generation is 



m(t) 



Pit) 
S(t) 



(2 - q)P 
l-(l-q)Pf 



(A3) 



Thus the ^-exponential function is obtained from a mortality 
function m(t) which is the inverse of a linear function of t. For 
fit «; 1, m(t) * (2 - g)/?, which is a small constant. For fit » 1, 
m{t) ss (2 - <20/(l - g)f, which is inversely proportional to t. 
This is reasonable in an ecological sense since species that have 
already existed a long time may be expected to also exist far into 
the future. 



Appendix B. Time series for lower mutation rates 

Although the mutation rates used in the body of the 
manuscript are high, the time series are intermittent even with 
lower mutation rates. Some examples of the time series with 
other parameter sets are shown in Fig. IB. 91 These time series 
are clearly different from the ones for the migration models. 
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